set.seed(5168)
## responsibility attribution ##

library(foreign)
library(mnormt)
library(ordinal)

data <- read.dta('Denmark (2014).dta')

matrix(c(1,2),1,2,byrow=TRUE)     #See what the matrix looks like
layout(matrix(c(1,2),1,2,byrow=TRUE), widths=c(1.2,1.2,1.2), heights=c(1.5,1.5,1.5,1.5,1.5), TRUE)  


######################
###GET COEFFICIENTS###
######################

model <- clmm(as.factor(Resp_attribution) ~ DK_V2_Prime_minister + DK_V2_Cabinet + DK_V2_Support + DK_V2_Opposition + DK_Perc_seats + DK_V2_Prime_minister * DK_Perc_seats + DK_V2_Cabinet * DK_Perc_seats + DK_V2_Support * DK_Perc_seats + DK_V2_Opposition*DK_Perc_seats + (1 | Party) + (1 | pid), data = data, HESS = TRUE, link = 'probit')
summary(model)

pars <- as.matrix(model$coefficients)
vce <- vcov(model)
vce <- vce[-c(14:15), ]       
vce <- vce[, -c(14:15)]       

simCoef <- rmnorm(1000, as.vector(pars),as.matrix(vce))

###############################################################################
### Marginal effect of perceived support party relative to perceived partner###
###############################################################################

pmHi <- c()
pmLo <- c()
pmMed <- c()

partHi <- c()
partLo <- c()
partMed <- c()

supHi <- c()
supLo <- c()
supMed <- c()

oppHi <- c()
oppLo <- c()
oppMed <- c()

otherHi <- c()
otherLo <- c()
otherMed <- c()

part_pred <- c()
sup_pred <- c()

marg_effect <- c()
margeff_Hi <- c()
margeff_Lo <- c()
margeff_Med <- c()



for(i in 0:20){
  ## Support party
  
  agg <- simCoef[, 5] * 0 + 
    simCoef[, 6] * 0 + 
    simCoef[, 7] * 1 + 
    simCoef[, 8] * 0 + 
    simCoef[, 9] * i + 
    simCoef[, 10] * 0 + 
    simCoef[, 11] * 0 + 
    simCoef[, 12] * i + 
    simCoef[, 13] * 0 
  
  pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
  pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
  pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
  pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
  pred5 <- 1 - pred1 - pred2 - pred3 - pred4
  
  pred <- pred5 * 5 + 
    pred4 * 4 + 
    pred3 * 3 + 
    pred2 * 2 + 
    pred1
  
  supHi[i + 1] <- quantile(pred, 0.975)
  supLo[i + 1] <- quantile(pred, 0.025)
  supMed[i + 1] <- quantile(pred, 0.5)
  
  sup_pred <- pred5 * 5 + 
    pred4 * 4 + 
    pred3 * 3 + 
    pred2 * 2 + 
    pred1
  
  ## partner party
  agg <- simCoef[, 5] * 0 + 
    simCoef[, 6] * 1 + 
    simCoef[, 7] * 0 + 
    simCoef[, 8] * 0 + 
    simCoef[, 9] * i + 
    simCoef[, 10] * 0 + 
    simCoef[, 11] * i + 
    simCoef[, 12] * 0 + 
    simCoef[, 13] * 0 
  
  pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
  pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
  pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
  pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
  pred5 <- 1 - pred1 - pred2 - pred3 - pred4
  
  pred <- pred5 * 5 + 
    pred4 * 4 + 
    pred3 * 3 + 
    pred2 * 2 + 
    pred1
  
  partHi[i + 1] <- quantile(pred, 0.975)
  partLo[i + 1] <- quantile(pred, 0.025)
  partMed[i + 1] <- quantile(pred, 0.5)
  
  part_pred <- pred5 * 5 + 
    pred4 * 4 + 
    pred3 * 3 + 
    pred2 * 2 + 
    pred1
  
  
  #get marginal effect 
  marg_effect = sup_pred - part_pred
  margeff_Hi[i + 1] = quantile(marg_effect, 0.975)
  margeff_Lo[i + 1] = quantile(marg_effect, 0.025)
  margeff_Med[i + 1] = quantile(marg_effect, 0.5)
}

## now paint a pretty picture.... ##

#seats sequence
seats <- seq(0, 20, 1)

#Now create the figure
par(mar=c(1,4.5,1.5,1.5))      
plot(1, 1, type = 'n',
     xlim = c(0, 20),
     ylim = c(-1, 1),
     xlab = '',
     ylab = '',
     main = 'Marginal effect of support role \nrelative to partner role', 
     cex.main=0.8,
     axes = FALSE)
	axis(1,at=seq(0,20,10),labels=T)         
axis(2)
polygon(c(seats, rev(seats)), c(margeff_Hi, rev(margeff_Lo)), col = rgb(red = 0.1, 0.1, 0, 0.2), border = FALSE)
lines(seats, margeff_Med, col = rgb(0.1, 0.1, 0, 1), lwd = 3)
abline(h=0,col="red")

mtext("Change in \nresponsibility attribution", side=2, line=2.1, cex=1)    
mtext("Perceived seat %", side=1, line=2.1, cex=1)         

##################################################################################
### Marginal effect of perceived support party relative to perceived opposition###
##################################################################################

pmHi <- c()
pmLo <- c()
pmMed <- c()

partHi <- c()
partLo <- c()
partMed <- c()

supHi <- c()
supLo <- c()
supMed <- c()

oppHi <- c()
oppLo <- c()
oppMed <- c()

otherHi <- c()
otherLo <- c()
otherMed <- c()

part_pred <- c()
sup_pred <- c()

marg_effect <- c()
margeff_Hi <- c()
margeff_Lo <- c()
margeff_Med <- c()



for(i in 0:20){
  ## Support party 
  
  agg <- simCoef[, 5] * 0 + 
    simCoef[, 6] * 0 + 
    simCoef[, 7] * 1 + 
    simCoef[, 8] * 0 + 
    simCoef[, 9] * i + 
    simCoef[, 10] * 0 + 
    simCoef[, 11] * 0 + 
    simCoef[, 12] * i + 
    simCoef[, 13] * 0 
  
  pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
  pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
  pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
  pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
  pred5 <- 1 - pred1 - pred2 - pred3 - pred4
  
  pred <- pred5 * 5 + 
    pred4 * 4 + 
    pred3 * 3 + 
    pred2 * 2 + 
    pred1
  
  supHi[i + 1] <- quantile(pred, 0.975)
  supLo[i + 1] <- quantile(pred, 0.025)
  supMed[i + 1] <- quantile(pred, 0.5)
  
  sup_pred <- pred5 * 5 + 
    pred4 * 4 + 
    pred3 * 3 + 
    pred2 * 2 + 
    pred1
  
  ## Opposition parties
  agg <- simCoef[, 5] * 0 + 
    simCoef[, 6] * 0 + 
    simCoef[, 7] * 0 + 
    simCoef[, 8] * 1 + 
    simCoef[, 9] * i + 
    simCoef[, 10] * 0 + 
    simCoef[, 11] * 0 + 
    simCoef[, 12] * 0 + 
    simCoef[, 13] * i 
  
  pred1 <- pnorm(simCoef[, 1], agg, sd = 1)
  pred2 <- pnorm(simCoef[, 2], agg, sd = 1) - pred1
  pred3 <- pnorm(simCoef[, 3], agg, sd = 1) - pred2 - pred1
  pred4 <- pnorm(simCoef[, 4], agg, sd = 1) - pred3 - pred2 - pred1
  pred5 <- 1 - pred1 - pred2 - pred3 - pred4
  
  pred <- pred5 * 5 + 
    pred4 * 4 + 
    pred3 * 3 + 
    pred2 * 2 + 
    pred1
  
  partHi[i + 1] <- quantile(pred, 0.975)
  partLo[i + 1] <- quantile(pred, 0.025)
  partMed[i + 1] <- quantile(pred, 0.5)
  
  part_pred <- pred5 * 5 + 
    pred4 * 4 + 
    pred3 * 3 + 
    pred2 * 2 + 
    pred1
  
  
  #get marginal effect 
  marg_effect = sup_pred - part_pred
  margeff_Hi[i + 1] = quantile(marg_effect, 0.975)
  margeff_Lo[i + 1] = quantile(marg_effect, 0.025)
  margeff_Med[i + 1] = quantile(marg_effect, 0.5)
}

## now paint a pretty picture.... ##

#seats sequence
seats <- seq(0, 20, 1)

#Now create the figure
par(mar=c(1,3,1.5,2))      
plot(1, 1, type = 'n',
     xlim = c(0, 20),
     ylim = c(-1, 1),
     xlab = '',
     ylab = '',
     main = 'Marginal effect of support role \nrelative to opposition role', 
     cex.main=0.8,
     axes = FALSE)
axis(1,at=seq(0,20,10),labels=T)          
axis(2)
polygon(c(seats, rev(seats)), c(margeff_Hi, rev(margeff_Lo)), col = rgb(red = 0.1, 0.1, 0, 0.2), border = FALSE)
lines(seats, margeff_Med, col = rgb(0.1, 0.1, 0, 1), lwd = 3)
abline(h=0,col="red")

mtext("Perceived seat %", side=1, line=2.1, cex=1)        

dev.copy(png,'Figure 6b DK 2014 - marginal effects of party types over seat share.png', height = 439, width = 600)
dev.off()
